This logbook will answer the Exercises of the module “Tree phenology
analysis in R”. It will start with chapter 3, ends with chapter 35 and
will only shortly tell if a chapter is left out (like chapter 17).
The logbook is split into the following 3 parts, Climate Scenarios
(chapter 3-19), PLS regression (chapter 20-28) and the phenoFlex model
(chapter 29-35).
This is part 1 about climate scenarios.
Put yourself in the place of a breeder who wants to calculate the
temperature requirements of a newly released cultivar. Which method will
you use to calculate the chilling and forcing periods?
Which are the advantages (name 2) of the BBCH scale compared with early scales?
Growth stages are easily recognizable under field conditions.
Growth stages are graded in the order of appearance.
Classify the following phenological stages of sweet cherry according to the BBCH scale.
The pictures are described from left to right. All of the classifications are based on the given pictures and I assume that the image section shown corresponds to the most represented of the tree.
Picture 01: Phenological stages
of a cherry tree
The reason are mainly human activities like deforestation, fossil fuel combustion and industrial processes
Greenhouse gases trap heat in the atmosphere
Can absorb only radiation of certain wavelengths
Short-wave radiation from Earth surface is absorbed
Recently the mean temperatures are globally rising clearly and very drastically
Especially in the last 50 to 100 years
This an other factors have a great influence in our climate and weather (regional and global)
We will and already do experience more extreme weather (like more dry periods, extreme rain, especially mild and extreme winter)
Most important is the rapid warming (as conclusion)
RCP -> Representative Concentration Pathways
RCPs project future greenhouse gas concentrations (NOT emissions) based on different climate change scenarios until the year 2100
The different scenarios are labeled after the possible radioative forcing values in the year 2100
Identifies regions where current conditions resemble projected future climates to predict potential species shifts.
Does not consider non-climatic factors (e.g. land use, human impact) that influence species survival.
Climate models often introduce biases that must be corrected before using their projections.
Different greenhouse gas emission scenarios (e.g. RCP2.6, RCP4.5, RCP8.5) result in a range of possible future conditions, adding uncertainty.
Statistical methods, such as ensemble modeling and probabilistic approaches, help quantify and communicate these uncertainties.
Data collection and Preprocessing
Chill Model Application
Climate Scenario Analysis
Validation and Uncertainty Assessment
Compare modeled projections with historical observations to evaluate accuracy.
Apply statistical techniques to quantify uncertainties and assess model reliability.
Use ensemble approaches to minimize errors and account for variability between models.
Visualization and Interpretation
Present results using graphs, maps, and reports for easy interpretation.
Provide insights for stakeholders, such as farmers and policymakers, to support decision-making.
Communicate potential risks and adaptation strategies based on projected changes.
Write a basic function that calculates warm hours (>25° C).
Basic structure:
# read in Dataset (e.g. Winters_hours_gaps)
dataset[ , new_column_name] <- dataset$temperature_column > 25
Apply this function to the Winters_hours_gaps dataset.
Winters_hours_gaps <- Winters_hours_gaps
hourtemps <- Winters_hours_gaps[ , c("Year", "Month", "Day", "Hour", "Temp")]
hourtemps[ , "warm_hour"] <- hourtemps$Temp > 25
hourtemps[2503:2507, ]
## Year Month Day Hour Temp warm_hour
## 2503 2008 6 15 16 27.604 TRUE
## 2504 2008 6 15 17 26.989 TRUE
## 2505 2008 6 15 18 24.436 FALSE
## 2506 2008 6 15 19 23.088 FALSE
## 2507 2008 6 15 20 19.555 FALSE
Start_Date <- which(hourtemps$Year == 2008 &
hourtemps$Month == 6 &
hourtemps$Day == 1 &
hourtemps$Hour == 12)
End_date <- which(hourtemps$Year == 2008 &
hourtemps$Month == 6 &
hourtemps$Day == 30 &
hourtemps$Hour ==12)
sum(hourtemps$warm_hour[Start_Date:End_date])
## [1] 250
# The result is 250 and will be the control for the following function with the same start and end
sum_warmH <- function(WHourly,
startYear,
startMonth,
startDay,
startHour,
endYear,
endMonth,
endDay,
endHour)
{WHourly[,"warm_hour2"] <- WHourly$Temp > 25
Start_row <- which(hourtemps$Year == startYear &
hourtemps$Month == startMonth &
hourtemps$Day == startDay &
hourtemps$Hour == startHour)
End_row <- which(hourtemps$Year == endYear &
hourtemps$Month == endMonth &
hourtemps$Day == endDay &
hourtemps$Hour == endHour)
WHs <- sum(WHourly$warm_hour[Start_row:End_row])
return(WHs)}
sum_warmH(hourtemps,
startYear = 2008,
startMonth = 6,
startDay = 1,
startHour = 12,
endYear = 2008,
endMonth = 6,
endDay = 30,
endHour = 12)
## [1] 250
hourtemps <- Winters_hours_gaps[ , c("Year", "Month", "Day", "Hour", "Temp")]
chilling_l7 <- chilling(make_JDay(hourtemps),
Start_JDay = 90,
End_JDay = 100)
chilling_l7
## Season End_year Season_days Data_days Perc_complete Chilling_Hours
## 1 2007/2008 2008 11 11 100 40
## Utah_Model Chill_portions GDH
## 1 15.5 2.009147 2406.52
df_l7 = data.frame(
lower = c(-1000, 1, 2, 3, 4, 5, 6),
upper = c( 1, 2, 3, 4, 5, 6, 1000),
weight = c( 0, 1, 2, 3, 2, 1, 0))
step_model_l7 <- step_model(HourTemp=hourtemps$Temp, df=df_l7)
TResp_l7 <- tempResponse(make_JDay(hourtemps),
Start_JDay = 90,
End_JDay = 100,
models = list(Chill_Portions = Dynamic_Model,
GDH = GDH))
TResp_l7
## Season End_year Season_days Data_days Perc_complete Chill_Portions GDH
## 1 2007/2008 2008 11 11 100 2.009147 2406.52
Days_l8 <- daylength(latitude = 51.4, JDay = 1:365)
Days_df_l8 <- data.frame(JDay = 1:365,
Sunrise = Days_l8$Sunrise,
Sunset = Days_l8$Sunset,
Daylength = Days_l8$Daylength)
Days_df_l8 <- pivot_longer(Days_df_l8, cols=c(Sunrise:Daylength))
ggplot(Days_df_l8, aes(JDay, value)) +
geom_line(lwd = 1.5) +
facet_grid(cols = vars(name)) +
ylab("Time of Day / Daylength [hours]") +
xlab("Day of the year [JDay]") +
theme_bw(base_size = 20)
hourtemp_KA <- stack_hourly_temps(KA_weather, latitude=50.4)
empi_curve_l8 <- Empirical_daily_temperature_curve(Winters_hours_gaps)
ggplot(data = empi_curve_l8[1:96, ], aes(Hour, Prediction_coefficient)) +
geom_line(lwd = 1.3, col = "red") +
facet_grid(rows = vars(Month)) +
xlab("Hour of the day") +
ylab("Prediction coefficient") +
theme_bw(base_size = 20)
coeffs_l8 <- Empirical_daily_temperature_curve(Winters_hours_gaps)
Winters_daily_l8 <- make_all_day_table(Winters_hours_gaps, input_timestep = "hour")
Winter_hours_l8 <- Empirical_hourly_temperatures(Winters_daily_l8, coeffs_l8)
# simplify the data for easier use
Winter_hours_l8 <- Winter_hours_l8[ , c("Year", "Month", "Day", "Hour", "Temp")]
colnames(Winter_hours_l8)[ncol(Winter_hours_l8)] <- "Temp_empirical"
Winters_ideal_l8 <- stack_hourly_temps(Winters_daily_l8, latitude = 38.5)$hourtemps
Winters_ideal_l8 <- Winters_ideal_l8[ , c("Year", "Month", "Day", "Hour", "Temp")]
colnames(Winters_ideal_l8)[ncol(Winters_ideal_l8)] <- "Temp_ideal"
# create the "triangular" data set to compare in the end
Winters_triangle_l8 <- Winters_daily_l8
Winters_triangle_l8[ , "Hour"] <- 0
Winters_triangle_l8$Hour[nrow(Winters_triangle_l8)] <- 23
Winters_triangle_l8[ , "Temp"] <- 0
Winters_triangle_l8 <- make_all_day_table(Winters_triangle_l8, timestep = "hour")
colnames(Winters_triangle_l8)[ncol(Winters_triangle_l8)] <- "Temp_triangular"
# the following loop will fill in the daily Tmin and Tmax values for every hour
for (i in 2:nrow(Winters_triangle_l8))
{if (is.na(Winters_triangle_l8$Tmin[i]))
Winters_triangle_l8$Tmin[i] <- Winters_triangle_l8$Tmin[i - 1]
if (is.na(Winters_triangle_l8$Tmax[i]))
Winters_triangle_l8$Tmax[i] <- Winters_triangle_l8$Tmax[i - 1]}
Winters_triangle_l8$Temp_triangular <- NA
# Tmin is fixated to the 6th hour of each day
Winters_triangle_l8$Temp_triangular[which(Winters_triangle_l8$Hour == 6)] <-
Winters_triangle_l8$Tmin[which(Winters_triangle_l8$Hour == 6)]
# Tmax is fixated to the 18th hour of each day
Winters_triangle_l8$Temp_triangular[which(Winters_triangle_l8$Hour == 18)] <-
Winters_triangle_l8$Tmax[which(Winters_triangle_l8$Hour == 18)]
# the following loop will fill in the gaps between the min and max data in a linear way
Winters_triangle_l8$Temp_triangular <- interpolate_gaps(Winters_triangle_l8$Temp_triangular)$interp
Winters_triangle_l8 <- Winters_triangle_l8[ , c("Year", "Month", "Day", "Hour", "Temp_triangular")]
# merge all created data frames for easy display
Winters_temps_l8 <- merge(Winters_hours_gaps, Winter_hours_l8,
by = c("Year", "Month", "Day", "Hour"))
Winters_temps_l8 <- merge(Winters_temps_l8, Winters_triangle_l8,
by = c("Year", "Month", "Day", "Hour"))
Winters_temps_l8 <- merge(Winters_temps_l8, Winters_ideal_l8,
by = c("Year", "Month", "Day", "Hour"))
# convert the date into R's date format and reorganize the data frame
Winters_temps_l8[, "DATE"] <-
ISOdate(Winters_temps_l8$Year,
Winters_temps_l8$Month,
Winters_temps_l8$Day,
Winters_temps_l8$Hour)
Winters_temps_to_plot_l8 <-
Winters_temps_l8[, c("DATE",
"Temp",
"Temp_empirical",
"Temp_triangular",
"Temp_ideal")]
Winters_temps_to_plot_l8 <- Winters_temps_to_plot_l8[100:200, ]
Winters_temps_to_plot_l8 <- pivot_longer(Winters_temps_to_plot_l8, cols=Temp:Temp_ideal)
colnames(Winters_temps_to_plot_l8) <- c("DATE", "Method", "Temperature")
ggplot(data = Winters_temps_to_plot_l8, aes(DATE, Temperature, colour = Method)) +
geom_line(lwd = 1.3) + ylab("Temperature (°C)") + xlab("Date")
# Convert the data set into a tibble
hourtemps <- Winters_hours_gaps
hourtemps %>% as_tibble() %>% summary()
## Year Month Day Hour
## Min. :2008 Min. : 3.000 Min. : 1.00 Min. : 0.0
## 1st Qu.:2008 1st Qu.: 5.000 1st Qu.: 8.00 1st Qu.: 6.0
## Median :2008 Median : 7.000 Median :15.00 Median :11.0
## Mean :2008 Mean : 6.722 Mean :15.53 Mean :11.5
## 3rd Qu.:2008 3rd Qu.: 9.000 3rd Qu.:23.00 3rd Qu.:17.0
## Max. :2008 Max. :11.000 Max. :31.00 Max. :23.0
##
## Temp_gaps Temp
## Min. : 1.967 Min. : 1.967
## 1st Qu.:13.016 1st Qu.:13.257
## Median :17.915 Median :18.010
## Mean :18.419 Mean :18.644
## 3rd Qu.:23.569 3rd Qu.:23.857
## Max. :38.170 Max. :38.365
## NA's :3374
# Select only the top 10 rows of the data set
hourtemps_l9 <- as_tibble(hourtemps[1:10, ])
# Convert the tibble to a long format, with separate rows for Temp_gaps and Temp
hourtemps_long_l9 <- hourtemps_l9 %>% pivot_longer(cols = Temp:Temp_gaps)
# Use ggplot2 to plot Temp_gaps and temp as facets (point or line plot)
ggplot(hourtemps_l9, aes(Hour, Temp_gaps)) + geom_point()
# Convert the data set back to the wide format
hourtemps_wide_l9 <- hourtemps_long_l9 %>% pivot_wider(names_from = name)
# Select only the following columns: Year, Month, Day and Temp
hourtemps_wide_l9 %>% select(c(Month, Day, Temp))
## # A tibble: 10 × 3
## Month Day Temp
## <int> <int> <dbl>
## 1 3 3 15.1
## 2 3 3 17.2
## 3 3 3 18.7
## 4 3 3 18.7
## 5 3 3 18.8
## 6 3 3 19.5
## 7 3 3 19.3
## 8 3 3 17.7
## 9 3 3 15.4
## 10 3 3 12.7
# Sort the data set by the Temp column, in descending order
hourtemps_wide_l9 %>% arrange(desc(Temp))
## # A tibble: 10 × 6
## Year Month Day Hour Temp Temp_gaps
## <int> <int> <int> <int> <dbl> <dbl>
## 1 2008 3 3 15 19.5 19.5
## 2 2008 3 3 16 19.3 19.3
## 3 2008 3 3 14 18.8 18.8
## 4 2008 3 3 12 18.7 18.7
## 5 2008 3 3 13 18.7 18.7
## 6 2008 3 3 17 17.7 17.7
## 7 2008 3 3 11 17.2 17.2
## 8 2008 3 3 18 15.4 15.4
## 9 2008 3 3 10 15.1 15.1
## 10 2008 3 3 19 12.7 12.7
HT_l9.2 <- Winters_hours_gaps
HT_l9.2$Temp_F <- NA
for (i in 1:nrow(HT_l9.2))
{HT_l9.2$Temp_F[i] <- (HT_l9.2$Temp[i] * 9/5) + 32}
head(HT_l9.2)
## Year Month Day Hour Temp_gaps Temp Temp_F
## 1 2008 3 3 10 15.127 15.127 59.2286
## 2 2008 3 3 11 17.153 17.153 62.8754
## 3 2008 3 3 12 18.699 18.699 65.6582
## 4 2008 3 3 13 18.699 18.699 65.6582
## 5 2008 3 3 14 18.842 18.842 65.9156
## 6 2008 3 3 15 19.508 19.508 67.1144
HT_l9.3 <- hourtemps
func_Far_l9.3 <- function(x)
{(x * 9/5) + 32}
head(HT_l9.3)
## Year Month Day Hour Temp_gaps Temp
## 1 2008 3 3 10 15.127 15.127
## 2 2008 3 3 11 17.153 17.153
## 3 2008 3 3 12 18.699 18.699
## 4 2008 3 3 13 18.699 18.699
## 5 2008 3 3 14 18.842 18.842
## 6 2008 3 3 15 19.508 19.508
for (i in 1:nrow(hourtemps))
{HT_l9.3$Temp_F[i] <- sapply(HT_l9.3$Temp[i], func_Far_l9.3)}
head(HT_l9.3)
## Year Month Day Hour Temp_gaps Temp Temp_F
## 1 2008 3 3 10 15.127 15.127 59.2286
## 2 2008 3 3 11 17.153 17.153 62.8754
## 3 2008 3 3 12 18.699 18.699 65.6582
## 4 2008 3 3 13 18.699 18.699 65.6582
## 5 2008 3 3 14 18.842 18.842 65.9156
## 6 2008 3 3 15 19.508 19.508 67.1144
HT_l9.2 <- hourtemps %>% mutate(Temp_F = (Temp * 9/5) + 32)
head(HT_l9.2)
## Year Month Day Hour Temp_gaps Temp Temp_F
## 1 2008 3 3 10 15.127 15.127 59.2286
## 2 2008 3 3 11 17.153 17.153 62.8754
## 3 2008 3 3 12 18.699 18.699 65.6582
## 4 2008 3 3 13 18.699 18.699 65.6582
## 5 2008 3 3 14 18.842 18.842 65.9156
## 6 2008 3 3 15 19.508 19.508 67.1144
station_list <- handle_gsod(action = "list_stations",
location = c(6.29, 51.40),
time_interval = c(1990, 2020))
station_list <- slice(station_list, 1:25)
head(station_list)
## # A tibble: 6 × 10
## chillR_code STATION.NAME CTRY Lat Long BEGIN END Distance
## <chr> <chr> <chr> <dbl> <dbl> <int> <int> <dbl>
## 1 06391099999 ARCEN AWS "NL" 51.5 6.2 19940901 20241127 12.8
## 2 10403099999 MOENCHENGLADBACH "GM" 51.2 6.5 19381001 19421031 23.6
## 3 10437499999 MONCHENGLADBACH "GM" 51.2 6.50 19960715 20241127 24.1
## 4 10405099999 LAARBRUCH (RAF) "GM" 51.6 6.15 19730102 19971231 24.3
## 5 10000199999 NIEDERRHEIN "GM" 51.6 6.14 20050101 20241127 24.7
## 6 99860099999 ARRC MOENCHENGLADBA "" 51.2 6.21 20011015 20020306 24.7
## # ℹ 2 more variables: Overlap_years <dbl>, Perc_interval_covered <dbl>
weather_WES <- handle_gsod(action = "download_weather",
location = station_list$chillR_code[1],
time_interval = c(1990, 2020))
cleaned_weather_WES <- handle_gsod(weather_WES)
cleaned_weather_WES[[1]][1000:1010,]
write.csv(station_list, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/station_list_Wesel.csv", row.names = FALSE)
write.csv(weather_WES[[1]], "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_raw_weather.csv", row.names = FALSE)
write.csv(cleaned_weather_WES[[1]], "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_chillR_weather.csv", row.names = FALSE)
Before the answers of the exercises for this chapter start, there will be dates added to the dataset because the first 4 years are missing for a reason. So the get added here manually to answer all following exercises correctly.
cleaned_weather_WES <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_chillR_weather.csv")
start_date <- as.Date("1990-01-01")
end_date <- as.Date("1993-12-31")
date_sequence <- seq.Date(start_date, end_date, by = "day")
new_rows <- data.frame(
Date = as.POSIXct(paste(date_sequence, "12:00:00")),
Year = as.numeric(format(date_sequence, "%Y")),
Month = as.numeric(format(date_sequence, "%m")),
Day = as.numeric(format(date_sequence, "%d")),
Tmin = NA,
Tmax = NA,
Tmean = NA,
Prec = NA
)
cleaned_weather_WES <- rbind(new_rows, cleaned_weather_WES)
Wesel_QC <- fix_weather(cleaned_weather_WES)$QC
table(is.na(cleaned_weather_WES$Tmin))
##
## FALSE TRUE
## 9022 2301
Wesel_QC
## Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1 1989/1990 1990 365 365 365 365
## 2 1990/1991 1991 365 365 365 365
## 3 1991/1992 1992 366 366 366 366
## 4 1992/1993 1993 365 365 365 365
## 5 1993/1994 1994 365 365 270 270
## 6 1994/1995 1995 365 365 161 161
## 7 1995/1996 1996 366 366 128 128
## 8 1996/1997 1997 365 365 152 152
## 9 1997/1998 1998 365 365 46 46
## 10 1998/1999 1999 365 365 6 6
## 11 1999/2000 2000 366 366 0 0
## 12 2000/2001 2001 365 365 31 31
## 13 2001/2002 2002 365 365 4 4
## 14 2002/2003 2003 365 365 2 2
## 15 2003/2004 2004 366 366 2 2
## 16 2004/2005 2005 365 365 11 11
## 17 2005/2006 2006 365 365 2 2
## 18 2006/2007 2007 365 365 4 4
## 19 2007/2008 2008 366 366 5 5
## 20 2008/2009 2009 365 365 0 0
## 21 2009/2010 2010 365 365 0 0
## 22 2010/2011 2011 365 365 0 0
## 23 2011/2012 2012 366 366 4 4
## 24 2012/2013 2013 365 365 5 5
## 25 2013/2014 2014 365 365 2 2
## 26 2014/2015 2015 365 365 2 2
## 27 2015/2016 2016 366 366 0 0
## 28 2016/2017 2017 365 365 0 0
## 29 2017/2018 2018 365 365 2 2
## 30 2018/2019 2019 365 365 0 0
## 31 2019/2020 2020 366 366 1 1
## Incomplete_days Perc_complete
## 1 365 0.0
## 2 365 0.0
## 3 366 0.0
## 4 365 0.0
## 5 270 26.0
## 6 161 55.9
## 7 128 65.0
## 8 152 58.4
## 9 46 87.4
## 10 6 98.4
## 11 0 100.0
## 12 31 91.5
## 13 4 98.9
## 14 2 99.5
## 15 2 99.5
## 16 11 97.0
## 17 2 99.5
## 18 4 98.9
## 19 5 98.6
## 20 0 100.0
## 21 0 100.0
## 22 0 100.0
## 23 4 98.9
## 24 5 98.6
## 25 2 99.5
## 26 2 99.5
## 27 0 100.0
## 28 0 100.0
## 29 2 99.5
## 30 0 100.0
## 31 1 99.7
station_list <- handle_gsod(action = "list_stations",
location = c(6.29, 51.40),
time_interval = c(1990, 2020))
station_list <- slice(station_list, 1:25)
station_list
## # A tibble: 25 × 10
## chillR_code STATION.NAME CTRY Lat Long BEGIN END Distance
## <chr> <chr> <chr> <dbl> <dbl> <int> <int> <dbl>
## 1 06391099999 ARCEN AWS "NL" 51.5 6.2 19940901 20241127 12.8
## 2 10403099999 MOENCHENGLADBACH "GM" 51.2 6.5 19381001 19421031 23.6
## 3 10437499999 MONCHENGLADBACH "GM" 51.2 6.50 19960715 20241127 24.1
## 4 10405099999 LAARBRUCH (RAF) "GM" 51.6 6.15 19730102 19971231 24.3
## 5 10000199999 NIEDERRHEIN "GM" 51.6 6.14 20050101 20241127 24.7
## 6 99860099999 ARRC MOENCHENGLADBA "" 51.2 6.21 20011015 20020306 24.7
## 7 10401099999 BRUGGEN (RAF) "GM" 51.2 6.13 19730102 19971231 24.8
## 8 06385099999 DE PEEL(NAFB) "NL" 51.6 5.93 19730402 19940611 29.9
## 9 10402099999 WILDENRATH(GAFB) "GM" 51.1 6.22 19730101 20030508 31.9
## 10 10404399999 KALKAR (MIL COMM) "GM" 51.7 6.17 19820927 19870304 32.6
## # ℹ 15 more rows
## # ℹ 2 more variables: Overlap_years <dbl>, Perc_interval_covered <dbl>
The stations 3, 5, 11 are suitable.
patch_weather <- handle_gsod(action = "download_weather",
location = as.character(station_list$chillR_code[c(3, 5, 11)]),
time_interval = c(1990,2020)) %>% handle_gsod()
## Loading data for 31 years from station 'NIEDERRHEIN'
## ================================================================================
##
## Loading data for 31 years from station 'DUSSELDORF'
## ================================================================================
##
## Loading data for 31 years from station 'MONCHENGLADBACH'
## ================================================================================
patched <- patch_daily_temperatures(weather = cleaned_weather_WES,
patch_weather = patch_weather,
max_mean_bias = 1,
max_stdev_bias = 2)
patched$statistics[[1]]
## mean_bias stdev_bias filled gaps_remain
## Tmin -0.703 1.151 36 2265
## Tmax 0.251 0.731 36 2265
patched$statistics[[2]]
## mean_bias stdev_bias filled gaps_remain
## Tmin -0.541 1.293 1098 1167
## Tmax 0.078 1.036 1098 1167
patched$statistics[[3]]
## mean_bias stdev_bias filled gaps_remain
## Tmin -1.206 1.705 0 1167
## Tmax 0.220 1.015 256 911
station_list
## # A tibble: 25 × 10
## chillR_code STATION.NAME CTRY Lat Long BEGIN END Distance
## <chr> <chr> <chr> <dbl> <dbl> <int> <int> <dbl>
## 1 06391099999 ARCEN AWS "NL" 51.5 6.2 19940901 20241127 12.8
## 2 10403099999 MOENCHENGLADBACH "GM" 51.2 6.5 19381001 19421031 23.6
## 3 10437499999 MONCHENGLADBACH "GM" 51.2 6.50 19960715 20241127 24.1
## 4 10405099999 LAARBRUCH (RAF) "GM" 51.6 6.15 19730102 19971231 24.3
## 5 10000199999 NIEDERRHEIN "GM" 51.6 6.14 20050101 20241127 24.7
## 6 99860099999 ARRC MOENCHENGLADBA "" 51.2 6.21 20011015 20020306 24.7
## 7 10401099999 BRUGGEN (RAF) "GM" 51.2 6.13 19730102 19971231 24.8
## 8 06385099999 DE PEEL(NAFB) "NL" 51.6 5.93 19730402 19940611 29.9
## 9 10402099999 WILDENRATH(GAFB) "GM" 51.1 6.22 19730101 20030508 31.9
## 10 10404399999 KALKAR (MIL COMM) "GM" 51.7 6.17 19820927 19870304 32.6
## # ℹ 15 more rows
## # ℹ 2 more variables: Overlap_years <dbl>, Perc_interval_covered <dbl>
patch_weather2 <- handle_gsod(action = "download_weather",
location = as.character(station_list$chillR_code[c(3, 4, 5, 7, 9, 11, 12, 13, 15, 16)]),
time_interval = c(1990,2020)) %>%
handle_gsod()
## Loading data for 31 years from station 'VOLKEL'
## ================================================================================
##
## Loading data for 31 years from station 'ELL AWS'
## ================================================================================
##
## Loading data for 31 years from station 'NIEDERRHEIN'
## ================================================================================
##
## Loading data for 31 years from station 'DUSSELDORF'
## ================================================================================
##
## Loading data for 31 years from station 'BRUGGEN (RAF)'
## ================================================================================
##
## Loading data for 31 years from station 'WILDENRATH(GAFB)'
## ================================================================================
##
## Loading data for 31 years from station 'KALKAR (MIL COMM)'
## ================================================================================
##
## Loading data for 31 years from station 'LAARBRUCH (RAF)'
## ================================================================================
##
## Loading data for 31 years from station 'ESSEN/MULHEIM'
## ================================================================================
##
## Loading data for 31 years from station 'MONCHENGLADBACH'
## ================================================================================
patched_2 <- patch_daily_temperatures(weather = cleaned_weather_WES,
patch_weather = patch_weather2,
max_mean_bias = 1,
max_stdev_bias = 2)
patched_2$statistics[[1]]
## mean_bias stdev_bias filled gaps_remain
## Tmin 0.665 1.179 2286 15
## Tmax 0.165 0.873 2286 15
patched_2$statistics[[6]]
## mean_bias stdev_bias filled gaps_remain
## Tmin 0.459 1.428 0 4
## Tmax -0.109 1.217 0 4
patched_2$statistics[[7]]
## mean_bias stdev_bias filled gaps_remain
## Tmin -0.302 1.258 3 1
## Tmax -0.018 1.211 3 1
post_patch_stats_2 <- fix_weather(patched_2)$QC
post_patch_stats_2
## Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1 1989/1990 1990 365 365 0 0
## 2 1990/1991 1991 365 365 0 0
## 3 1991/1992 1992 366 366 0 0
## 4 1992/1993 1993 365 365 0 0
## 5 1993/1994 1994 365 365 0 0
## 6 1994/1995 1995 365 365 0 0
## 7 1995/1996 1996 366 366 0 0
## 8 1996/1997 1997 365 365 0 0
## 9 1997/1998 1998 365 365 0 0
## 10 1998/1999 1999 365 365 1 1
## 11 1999/2000 2000 366 366 0 0
## 12 2000/2001 2001 365 365 0 0
## 13 2001/2002 2002 365 365 0 0
## 14 2002/2003 2003 365 365 0 0
## 15 2003/2004 2004 366 366 0 0
## 16 2004/2005 2005 365 365 0 0
## 17 2005/2006 2006 365 365 0 0
## 18 2006/2007 2007 365 365 0 0
## 19 2007/2008 2008 366 366 0 0
## 20 2008/2009 2009 365 365 0 0
## 21 2009/2010 2010 365 365 0 0
## 22 2010/2011 2011 365 365 0 0
## 23 2011/2012 2012 366 366 0 0
## 24 2012/2013 2013 365 365 0 0
## 25 2013/2014 2014 365 365 0 0
## 26 2014/2015 2015 365 365 0 0
## 27 2015/2016 2016 366 366 0 0
## 28 2016/2017 2017 365 365 0 0
## 29 2017/2018 2018 365 365 0 0
## 30 2018/2019 2019 365 365 0 0
## 31 2019/2020 2020 366 366 0 0
## Incomplete_days Perc_complete
## 1 0 100.0
## 2 0 100.0
## 3 0 100.0
## 4 0 100.0
## 5 0 100.0
## 6 0 100.0
## 7 0 100.0
## 8 0 100.0
## 9 0 100.0
## 10 1 99.7
## 11 0 100.0
## 12 0 100.0
## 13 0 100.0
## 14 0 100.0
## 15 0 100.0
## 16 0 100.0
## 17 0 100.0
## 18 0 100.0
## 19 0 100.0
## 20 0 100.0
## 21 0 100.0
## 22 0 100.0
## 23 0 100.0
## 24 0 100.0
## 25 0 100.0
## 26 0 100.0
## 27 0 100.0
## 28 0 100.0
## 29 0 100.0
## 30 0 100.0
## 31 0 100.0
patched_complete <- patched_2$weather %>% make_all_day_table()
Tmin_int <- interpolate_gaps(patched_complete[,"Tmin"])
patched_complete <- patched_complete %>% mutate(Tmin = Tmin_int$interp,
Tmin_interpolated = Tmin_int$missing)
Tmax_int <- interpolate_gaps(patched_complete[,"Tmax"])
patched_complete <- patched_complete %>% mutate(Tmax = Tmax_int$interp,
Tmax_interpolated = Tmax_int$missing)
fix_weather(patched_complete)$QC
## Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1 1989/1990 1990 365 365 0 0
## 2 1990/1991 1991 365 365 0 0
## 3 1991/1992 1992 366 366 0 0
## 4 1992/1993 1993 365 365 0 0
## 5 1993/1994 1994 365 365 0 0
## 6 1994/1995 1995 365 365 0 0
## 7 1995/1996 1996 366 366 0 0
## 8 1996/1997 1997 365 365 0 0
## 9 1997/1998 1998 365 365 0 0
## 10 1998/1999 1999 365 365 0 0
## 11 1999/2000 2000 366 366 0 0
## 12 2000/2001 2001 365 365 0 0
## 13 2001/2002 2002 365 365 0 0
## 14 2002/2003 2003 365 365 0 0
## 15 2003/2004 2004 366 366 0 0
## 16 2004/2005 2005 365 365 0 0
## 17 2005/2006 2006 365 365 0 0
## 18 2006/2007 2007 365 365 0 0
## 19 2007/2008 2008 366 366 0 0
## 20 2008/2009 2009 365 365 0 0
## 21 2009/2010 2010 365 365 0 0
## 22 2010/2011 2011 365 365 0 0
## 23 2011/2012 2012 366 366 0 0
## 24 2012/2013 2013 365 365 0 0
## 25 2013/2014 2014 365 365 0 0
## 26 2014/2015 2015 365 365 0 0
## 27 2015/2016 2016 366 366 0 0
## 28 2016/2017 2017 365 365 0 0
## 29 2017/2018 2018 365 365 0 0
## 30 2018/2019 2019 365 365 0 0
## 31 2019/2020 2020 366 366 0 0
## Incomplete_days Perc_complete
## 1 0 100
## 2 0 100
## 3 0 100
## 4 0 100
## 5 0 100
## 6 0 100
## 7 0 100
## 8 0 100
## 9 0 100
## 10 0 100
## 11 0 100
## 12 0 100
## 13 0 100
## 14 0 100
## 15 0 100
## 16 0 100
## 17 0 100
## 18 0 100
## 19 0 100
## 20 0 100
## 21 0 100
## 22 0 100
## 23 0 100
## 24 0 100
## 25 0 100
## 26 0 100
## 27 0 100
## 28 0 100
## 29 0 100
## 30 0 100
## 31 0 100
fix_QC <- fix_weather(patched_complete)$QC
fix_QC
## Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1 1989/1990 1990 365 365 0 0
## 2 1990/1991 1991 365 365 0 0
## 3 1991/1992 1992 366 366 0 0
## 4 1992/1993 1993 365 365 0 0
## 5 1993/1994 1994 365 365 0 0
## 6 1994/1995 1995 365 365 0 0
## 7 1995/1996 1996 366 366 0 0
## 8 1996/1997 1997 365 365 0 0
## 9 1997/1998 1998 365 365 0 0
## 10 1998/1999 1999 365 365 0 0
## 11 1999/2000 2000 366 366 0 0
## 12 2000/2001 2001 365 365 0 0
## 13 2001/2002 2002 365 365 0 0
## 14 2002/2003 2003 365 365 0 0
## 15 2003/2004 2004 366 366 0 0
## 16 2004/2005 2005 365 365 0 0
## 17 2005/2006 2006 365 365 0 0
## 18 2006/2007 2007 365 365 0 0
## 19 2007/2008 2008 366 366 0 0
## 20 2008/2009 2009 365 365 0 0
## 21 2009/2010 2010 365 365 0 0
## 22 2010/2011 2011 365 365 0 0
## 23 2011/2012 2012 366 366 0 0
## 24 2012/2013 2013 365 365 0 0
## 25 2013/2014 2014 365 365 0 0
## 26 2014/2015 2015 365 365 0 0
## 27 2015/2016 2016 366 366 0 0
## 28 2016/2017 2017 365 365 0 0
## 29 2017/2018 2018 365 365 0 0
## 30 2018/2019 2019 365 365 0 0
## 31 2019/2020 2020 366 366 0 0
## Incomplete_days Perc_complete
## 1 0 100
## 2 0 100
## 3 0 100
## 4 0 100
## 5 0 100
## 6 0 100
## 7 0 100
## 8 0 100
## 9 0 100
## 10 0 100
## 11 0 100
## 12 0 100
## 13 0 100
## 14 0 100
## 15 0 100
## 16 0 100
## 17 0 100
## 18 0 100
## 19 0 100
## 20 0 100
## 21 0 100
## 22 0 100
## 23 0 100
## 24 0 100
## 25 0 100
## 26 0 100
## 27 0 100
## 28 0 100
## 29 0 100
## 30 0 100
## 31 0 100
table(is.na(patched_complete$Tmin))
##
## FALSE
## 11323
write.csv(patched_complete, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv", row.names = FALSE)
patched_complete <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")
# we simulate here with the data form 1998 to 2009 100 times how the weather could have been in this time
Temp <- patched_complete %>%
temperature_generation(years = c(1998,2009),
sim_years = c(2001,2100))
## Warning: scenario doesn't contain named elements - consider using the following
## element names: 'data', 'reference_year','scenario_type','labels'
## Warning: setting scenario_type to 'relative'
## Warning: Reference year missing - can't check if relative temperature scenario
## is valid
Temperatures <- patched_complete %>% filter(Year %in% 1998:2009) %>%
cbind(Data_source = "observed") # cbind() adds more coloumns to the table
new_data <- Temp[[1]] %>%
select(Year, Month, Day, Tmin, Tmax) %>% # Nur die gewünschten Spalten auswählen
mutate(Data_source = "simulated") # Neue Spalte Data_source hinzufügen
# Daten an Temperatures anhängen
Temperatures <- bind_rows(Temperatures, new_data)
Temperatures <- Temperatures %>%
mutate(Date = as.Date(ISOdate(2000, Month, Day))) # mutate creates a coloumn from already existing data
ggplot(data = Temperatures,
aes(Date,
Tmin)) +
geom_smooth(aes(colour = factor(Year))) +
facet_wrap(vars(Data_source)) +
theme_bw(base_size = 20) +
theme(legend.position = "none") +
scale_x_date(date_labels = "%b")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(data = Temperatures,
aes(Date,
Tmax)) + # creates a empty coordinate system
geom_smooth(aes(colour = factor(Year))) + # creates a smoothing regression to make the data better obsevable (alternative is geom_line() then you getalso the nice colors but without the smoothing evect)
facet_wrap(vars(Data_source)) + #
theme_bw(base_size = 20) +
theme(legend.position = "none") + # deletes the legend (is here only confusing)
scale_x_date(date_labels = "%b") # shows on the xachsis only the month without the year
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
chill_observed <- Temperatures %>%
filter(Data_source == "observed") %>%
stack_hourly_temps(latitude = 50.4) %>%
chilling(Start_JDay = 305,
End_JDay = 59)
chill_simulated <- Temperatures %>%
filter(Data_source == "simulated") %>%
stack_hourly_temps(latitude = 51.4) %>%
chilling(Start_JDay = 305,
End_JDay = 59)
chill_comparison <-
cbind(chill_observed,
Data_source = "observed") %>%
rbind(cbind(chill_simulated,
Data_source = "simulated"))
chill_comparison_full_seasons <-
chill_comparison %>%
filter(Perc_complete == 100)
chill_comparison_full_seasons$Data_source <- factor(
chill_comparison_full_seasons$Data_source,
levels = c("simulated", "observed")
)
chill_comparison_full_seasons$Data_source <- factor(
chill_comparison_full_seasons$Data_source,
levels = c("simulated", "observed")
)
chill_comparison_full_seasons <- chill_comparison_full_seasons[order(
chill_comparison_full_seasons$Data_source), ]
ggplot(chill_comparison_full_seasons,
aes(x = Chill_portions)) +
geom_histogram(binwidth = 1, position = "identity",
aes(fill = factor(Data_source))) +
scale_fill_brewer(palette="Set2") +
theme_bw(base_size = 20) +
labs(fill = "Data source") +
xlab("Chill accumulation (Chill Portions)") +
ylab("Frequency")
chill_simulations <-
chill_comparison_full_seasons %>%
filter(Data_source == "simulated")
ggplot(chill_simulations,
aes(x = Chill_portions)) +
stat_ecdf(geom = "step",
lwd = 1.5,
col = "blue") +
ylab("Cumulative probability") +
xlab("Chill accumulation (in Chill Portions)") +
theme_bw(base_size = 20)
# filter only our observed data and calculate with the coordinates (latitude) the chilling day
chill_observed <- Temperatures %>%
filter(Data_source == "observed") %>%
stack_hourly_temps(latitude = 51.4) %>%
chilling(Start_JDay = 305,
End_JDay = 59)
### adding a freezing model - to make this for April, you'll also have to adjust the dates in the calculations
df <- data.frame(
lower= c(-1000, 0),
upper= c( 0, 1000),
weight=c( 1, 0))
freezing_hours <- function(x) step_model(x,df)
freezing_hours(c(1, 2, 4, 5, -10))
## [1] 0 0 0 0 1
chill_observed <- Temperatures %>%
filter(Data_source == "observed") %>%
stack_hourly_temps(latitude = 51.4) %>%
tempResponse(Start_JDay = 305, # more flexible than the chilling() function, we feed it our freezing function
End_JDay = 59,
models = list(Frost = freezing_hours,
Chill_portions = Dynamic_Model,
GDH = GDH))
####
# repeating the freezing function with the simulated data
chill_simulated <- Temperatures %>%
filter(Data_source == "simulated") %>%
stack_hourly_temps(latitude = 51.4) %>%
tempResponse(Start_JDay = 305,
End_JDay = 59,
models=list(Frost = freezing_hours,
Chill_portions = Dynamic_Model,
GDH = GDH))
# combine the freezing data for the observed and the simulated
chill_comparison <-
cbind(chill_observed,
Data_source = "observed") %>%
rbind(cbind(chill_simulated,
Data_source = "simulated"))
# sometimes there is data missing for the winter and the influence the data negative, so we remove them
chill_comparison_full_seasons <-
chill_comparison %>%
filter(Perc_complete == 100)
chill_comparison_full_seasons$Data_source <- factor(
chill_comparison_full_seasons$Data_source,
levels = c("simulated", "observed")
)
chill_comparison_full_seasons <- chill_comparison_full_seasons[order(
chill_comparison_full_seasons$Data_source), ]
# plot the data as a histogram
ggplot(chill_comparison_full_seasons,
aes(x = Frost)) +
geom_histogram(binwidth = 25, position = "identity",
aes(fill = factor(Data_source))) +
scale_fill_brewer(palette="Set2") +
theme_bw(base_size = 10) +
labs(fill = "Data source") +
xlab("Frost incidence during winter (hours)") +
ylab("Frequency")
chill_simulations <-
chill_comparison_full_seasons %>%
filter(Data_source == "simulated")
# shows how many frost hours are how likely to happen
ggplot(chill_simulations,
aes(x = Frost)) +
stat_ecdf(geom = "step", # accumulated density function
lwd = 1.5,
col = "blue") +
ylab("Cumulative probability") +
xlab("Frost incidence during winter (hours)") +
theme_bw(base_size = 20)
# Here's the amount of chill that is exceeded in 90% of all years.
quantile(chill_simulations$Chill_portions, 0.1)
## 10%
## 78.80063
# and here's the 50% confidence interval (25th to 75th percentile)
quantile(chill_simulations$Chill_portions, c(0.25, 0.75))
## 25% 75%
## 82.04394 86.97673
There are no tasks for this chapter.
patched_complete <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")
Temp <- patched_complete %>%
temperature_generation(years = c(1998,2005),
sim_years = c(2001,2100))
change_scenario <- data.frame(Tmin = rep(2,12),
Tmax = rep(2,12))
# change_scenario
Temp_2 <- temperature_generation(KA_weather,
years = c(1998,2005),
sim_years = c(2001,2100),
temperature_scenario = change_scenario)
Temperature_scenarios <- KA_weather %>%
filter(Year %in% 1998:2005) %>%
cbind(Data_source = "observed") %>%
rbind(Temp[[1]] %>%
select(c(Year, Month, Day, Tmin, Tmax)) %>%
cbind(Data_source = "simulated")
) %>%
rbind(Temp_2[[1]] %>%
select(c(Year, Month, Day, Tmin, Tmax)) %>%
cbind(Data_source = "Warming_2C")
) %>%
mutate(Date = as.Date(ISOdate(2000,
Month,
Day)))
ggplot(data = Temperature_scenarios,
aes(Date,Tmin)) +
geom_smooth(aes(colour = factor(Year))) +
facet_wrap(vars(Data_source)) +
theme_bw(base_size = 20) +
theme(legend.position = "none") +
scale_x_date(date_labels = "%b")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(data = Temperature_scenarios,
aes(Date,Tmax)) +
geom_smooth(aes(colour = factor(Year))) +
facet_wrap(vars(Data_source)) +
theme_bw(base_size = 20) +
theme(legend.position = "none") +
scale_x_date(date_labels = "%b")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
summary(patched_complete)
## YEARMODA DATE Date Year
## Min. :19900101 Length:11323 Length:11323 Min. :1990
## 1st Qu.:19971002 Class :character Class :character 1st Qu.:1997
## Median :20050702 Mode :character Mode :character Median :2005
## Mean :20050675 Mean :2005
## 3rd Qu.:20130402 3rd Qu.:2013
## Max. :20201231 Max. :2020
##
## Month Day Tmin Tmax
## Min. : 1.000 Min. : 1.00 Min. :-19.889 Min. :-10.278
## 1st Qu.: 4.000 1st Qu.: 8.00 1st Qu.: 2.276 1st Qu.: 9.222
## Median : 7.000 Median :16.00 Median : 6.722 Median : 15.111
## Mean : 6.523 Mean :15.73 Mean : 6.484 Mean : 15.112
## 3rd Qu.:10.000 3rd Qu.:23.00 3rd Qu.: 11.111 3rd Qu.: 21.000
## Max. :12.000 Max. :31.00 Max. : 23.278 Max. : 40.222
##
## Tmean Prec Tmin_source Tmax_source
## Min. :-15.222 Min. : 0.000 Length:11323 Length:11323
## 1st Qu.: 5.944 1st Qu.: 0.000 Class :character Class :character
## Median : 10.916 Median : 0.000 Mode :character Mode :character
## Mean : 10.784 Mean : 1.952
## 3rd Qu.: 15.944 3rd Qu.: 2.032
## Max. : 31.278 Max. :78.486
## NA's :2301 NA's :2301
## Tmin_interpolated Tmax_interpolated
## Mode :logical Mode :logical
## FALSE:11322 FALSE:11322
## TRUE :1 TRUE :1
##
##
##
##
scenario_1990 <- temperature_scenario_from_records(weather = patched_complete,
year = 1990)
temps_1990 <- temperature_generation(weather = patched_complete,
years = c(1990,2020),
sim_years = c(2001,2100),
temperature_scenario = scenario_1990)
scenario_2005 <- temperature_scenario_from_records(weather = patched_complete,
year = 2005)
relative_scenario <- temperature_scenario_baseline_adjustment(
baseline = scenario_2005,
temperature_scenario = scenario_1990)
temps_1990<-temperature_generation(weather = patched_complete,
years = c(1990,2020),
sim_years = c(2001,2100),
temperature_scenario = relative_scenario)
all_past_scenarios <- temperature_scenario_from_records(
weather = patched_complete,
year = c(1990,
1997,
2003,
2010))
adjusted_scenarios <- temperature_scenario_baseline_adjustment(
baseline = scenario_2005,
temperature_scenario = all_past_scenarios)
all_past_scenario_temps <- temperature_generation(
weather = patched_complete,
years = c(1990,2020),
sim_years = c(2001,2100),
temperature_scenario = adjusted_scenarios)
save_temperature_scenarios(all_past_scenario_temps, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_scenarios")
frost_model <- function(x)
step_model(x,
data.frame(
lower=c(-1000,0),
upper=c(0,1000),
weight=c(1,0)))
models <- list(Chill_Portions = Dynamic_Model,
GDH = GDH,
Frost_H = frost_model)
chill_hist_scenario_list <- tempResponse_daily_list(all_past_scenario_temps,
latitude = 51.4,
Start_JDay = 305,
End_JDay = 59,
models = models)
chill_hist_scenario_list <- lapply(chill_hist_scenario_list,
function(x) x %>%
filter(Perc_complete == 100))
save_temperature_scenarios(chill_hist_scenario_list, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_305_59")
scenarios <- names(chill_hist_scenario_list)[1:4]
all_scenarios <- chill_hist_scenario_list[[scenarios[1]]] %>%
mutate(scenario = as.numeric(scenarios[1]))
for (sc in scenarios[2:4])
all_scenarios <- all_scenarios %>%
rbind(chill_hist_scenario_list[[sc]] %>%
cbind(
scenario=as.numeric(sc))
) %>%
filter(Perc_complete == 100)
# Let's compute the actual 'observed' chill for comparison
actual_chill <- tempResponse_daily_list(patched_complete,
latitude=51.4,
Start_JDay = 305,
End_JDay = 59,
models)[[1]] %>%
filter(Perc_complete == 100)
ggplot(data = all_scenarios,
aes(scenario,
Chill_Portions,
fill = factor(scenario))) +
geom_violin() +
ylab("Chill accumulation (Chill Portions)") +
xlab("Scenario year") +
theme_bw(base_size = 15) +
ylim(c(0,90)) +
geom_point(data = actual_chill,
aes(End_year,
Chill_Portions,
fill = "blue"),
col = "blue",
show.legend = FALSE) +
scale_fill_discrete(name = "Scenario",
breaks = unique(all_scenarios$scenario))
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
write.csv(actual_chill,"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv", row.names = FALSE)
temperature_means <-
data.frame(Year = min(patched_complete$Year):max(patched_complete$Year),
Tmin = aggregate(patched_complete$Tmin,
FUN = "mean",
by = list(patched_complete$Year))[,2],
Tmax=aggregate(patched_complete$Tmax,
FUN = "mean",
by = list(patched_complete$Year))[,2]) %>%
mutate(runn_mean_Tmin = runn_mean(Tmin,15),
runn_mean_Tmax = runn_mean(Tmax,15))
Tmin_regression <- lm(Tmin~Year,
temperature_means)
Tmax_regression <- lm(Tmax~Year,
temperature_means)
temperature_means <- temperature_means %>%
mutate(regression_Tmin = Tmin_regression$coefficients[1]+
Tmin_regression$coefficients[2]*temperature_means$Year,
regression_Tmax = Tmax_regression$coefficients[1]+
Tmax_regression$coefficients[2]*temperature_means$Year
)
ggplot(temperature_means,
aes(Year,
Tmin)) +
geom_point() +
geom_line(data = temperature_means,
aes(Year,
runn_mean_Tmin),
lwd = 2,
col = "blue") +
geom_line(data = temperature_means,
aes(Year,
regression_Tmin),
lwd = 2,
col = "red") +
theme_bw(base_size = 15) +
ylab("Mean monthly minimum temperature (°C)")
ggplot(temperature_means,
aes(Year,
Tmax)) +
geom_point() +
geom_line(data = temperature_means,
aes(Year,
runn_mean_Tmax),
lwd = 2,
col = "blue") +
geom_line(data = temperature_means,
aes(Year,
regression_Tmax),
lwd = 2,
col = "red") +
theme_bw(base_size = 15) +
ylab("Mean monthly maximum temperature (°C)")
The main differences between RCPs (Representative Concentration Pathways) ans SSPs (Shared Socioeconomic Pathways) are the following:
Focus
RCPs are defined through greenhouse gas concentrations and radioactive forcing levels
SSPs include socioeconomic factors (e.g. population growth, energy use, policy decisions) along with emissions pathways
Drivers of Change
RCPs only consider atmospheric greenhouse gas concentrations and their effects on climate
SSPs include also social, economic and technological factors that influence emissions
Policy Relevance
RCPs provide no direct link to human actions or policy decisions
SSPs offer insights into how different societal choices impact climate change
Emissions Pathway Representation
RCPs directly define future greenhouse gas concentration levels
SSPs can be paired with radioactive forcing levels for more flexibility
RCPs are considered now as outdated and SSPs are now mainly used
location = c(6.29, 51.40)
area <- c(52.5, 5, 50.5, 7)
download_cmip6_ecmwfr(
scenarios = c("ssp126", "ssp245", "ssp370", "ssp585"),
area = area,
key = '9d522e37-4b90-46d1-8701-7deadd45032b',
model = 'default',
frequency = 'monthly',
variable = c('Tmin', 'Tmax'),
year_start = 2015,
year_end = 2100)
download_baseline_cmip6_ecmwfr(
area = area,
key = '9d522e37-4b90-46d1-8701-7deadd45032b',
model = 'match_downloaded',
frequency = 'monthly',
variable = c('Tmin', 'Tmax'),
year_start = 1986,
year_end = 2014,
month = 1:12)
library(ncdf4)
library(PCICt)
Wesel_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")
station <- data.frame(
station_name = c("ARCEN AWS"),
longitude = c(6.5),
latitude = c(51.8))
extracted <- extract_cmip6_data(stations = station)
## Unzipping files
## Extracting downloaded CMIP6 files
# head(extracted$`ssp126_AWI-CM-1-1-MR`)
change_scenarios <-
gen_rel_change_scenario(extracted,
scenarios = c(2050, 2085),
reference_period = c(1986:2014),
future_window_width = 30)
head(change_scenarios)
## # A tibble: 6 × 11
## location Month Tmax Tmin scenario start_year end_year scenario_year
## <chr> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 ARCEN AWS 1 1.90 1.95 ssp126 2035 2065 2050
## 2 ARCEN AWS 2 2.60 2.00 ssp126 2035 2065 2050
## 3 ARCEN AWS 3 1.40 1.04 ssp126 2035 2065 2050
## 4 ARCEN AWS 4 1.38 1.15 ssp126 2035 2065 2050
## 5 ARCEN AWS 5 1.87 1.56 ssp126 2035 2065 2050
## 6 ARCEN AWS 6 1.86 1.67 ssp126 2035 2065 2050
## # ℹ 3 more variables: reference_year <int>, scenario_type <chr>, labels <chr>
write.csv(change_scenarios, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/all_change_scenarios_ssp126.csv", row.names = FALSE)
change_scenarios <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/all_change_scenarios_ssp126.csv")
scen_list <- convert_scen_information(change_scenarios)
scen_frame <- convert_scen_information(scen_list)
# scen_list$"ARCEN AWS"$ssp126$`ACCESS-CM2`$'2050'
temps_2005 <- temperature_scenario_from_records(Wesel_temps,
2005)
temps_2000 <- temperature_scenario_from_records(Wesel_temps,
2000)
# temps_2005
# temps_2000
# Wesel_temps
base <- temperature_scenario_baseline_adjustment(temps_2005,
temps_2000)
scen_list <- convert_scen_information(change_scenarios,
give_structure = FALSE)
adjusted_list <- temperature_scenario_baseline_adjustment(base,
scen_list,
temperature_check_args =
list(scenario_check_thresholds = c(-5, 15)))
temps <- temperature_generation(Wesel_temps,
years = c(1990, 2020),
sim_years = c(2001, 2100),
adjusted_list,
temperature_check_args =
list( scenario_check_thresholds = c(-5, 15)))
save_temperature_scenarios(temps, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_future")
for(scen in 1:length(adjusted_list))
{
if(!file.exists(paste0("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_future_",
scen,"_",
names(adjusted_list)[scen],".csv")) )
{temp_temp <- temperature_generation(Wesel_temps,
years = c(1990, 2020),
sim_years = c(2001, 2100),
adjusted_list[scen],
temperature_check_args =
list( scenario_check_thresholds = c(-5, 15)))
write.csv(temp_temp[[1]],paste0("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_future_",scen,"_",names(adjusted_list)[scen],".csv"),
row.names=FALSE)
print(paste("Processed object",scen,"of", length(adjusted_list)))
}
}
temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_future_")
frost_model <- function(x)
step_model(x,
data.frame(
lower = c(-1000, 0),
upper = c(0, 1000),
weight = c(1, 0)))
models <- list(Chill_Portions = Dynamic_Model,
GDH = GDH,
Frost_H = frost_model)
chill_future_scenario_list <-
tempResponse_daily_list(temps,
latitude = 51.4,
Start_JDay = 305,
End_JDay = 59,
models = models)
chill_future_scenario_list <-
lapply(chill_future_scenario_list,
function(x) x %>%
filter(Perc_complete == 100))
save_temperature_scenarios(chill_future_scenario_list,
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
"Wesel_futurechill")
chill_future_scenario_list_frost <-
lapply(chill_future_scenario_list_frost,
function(x) x %>%
filter(Perc_complete == 100))
save_temperature_scenarios(chill_future_scenario_list_frost,
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
"Wesel_futurefrost")
Wesel_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")
observed_chill <- tempResponse_daily_list(Wesel_temps,
latitude = 51.4,
Start_JDay = 305,
End_JDay = 59,
models = models)
observed_chill <- lapply(observed_chill,
function(x) x %>%
filter(Perc_complete == 100))
write.csv(observed_chill, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv")
hist_temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_")
chill_hist_scenario_list <- tempResponse_daily_list(hist_temps,
latitude = 51.4,
Start_JDay = 305,
End_JDay = 59,
models = models)
chill_hist_scenario_list <- lapply(chill_hist_scenario_list,
function(x) x %>%
filter(Perc_complete == 100))
save_temperature_scenarios(chill_hist_scenario_list,
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
"Wesel_hist_chill_305_59")
chill_future_scenario_list <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_futurechill")
chill_hist_scenario_list<-load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_305_59")
observed_chill <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv")
chills <- make_climate_scenario(
metric_summary = chill_hist_scenario_list,
caption = "Historical",
historic_data = observed_chill,
time_series = TRUE)
plot_climate_scenarios(
climate_scenario_list = chills,
metric = "Chill_Portions",
metric_label = "Chill (Chill Portions)")
## [[1]]
## [1] "time series labels"
list_ssp <-
strsplit(names(chill_future_scenario_list),
'\\.') %>%
map(2) %>%
unlist()
list_gcm <-
strsplit(names(chill_future_scenario_list), '\\.') %>%
map(3) %>%
unlist()
list_time <-
strsplit(names(chill_future_scenario_list), '\\.') %>%
map(4) %>%
unlist()
SSPs <- c("ssp126", "ssp245", "ssp370", "ssp585")
Times <- c(2050, 2085)
for(SSP in SSPs)
for(Time in Times)
{
# find all scenarios for the ssp and time
chill <- chill_future_scenario_list[list_ssp == SSP &
list_time == Time]
names(chill) <- list_gcm[list_ssp == SSP &
list_time == Time]
if(SSP == "ssp126") SSPcaption <- "SSP1"
if(SSP == "ssp245") SSPcaption <- "SSP2"
if(SSP == "ssp370") SSPcaption <- "SSP3"
if(SSP == "ssp585") SSPcaption <- "SSP5"
chills <- chill %>%
make_climate_scenario(
caption = c(SSPcaption,
Time),
add_to = chills)
}
info_chill <-
plot_climate_scenarios(
climate_scenario_list = chills,
metric = "Chill_Portions",
metric_label = "Chill (Chill Portions)",
texcex = 1)
info_heat <-
plot_climate_scenarios(
climate_scenario_list = chills,
metric = "GDH",
metric_label = "Heat (Growing Degree Hours)",
texcex = 1)
info_frost <-
plot_climate_scenarios(
climate_scenario_list=chills,
metric="Frost_H",
metric_label="Frost incidence (hours)",
texcex=1)
info_chill[[2]]
## code Label
## 1 1 INM-CM4-8
## 2 2 MPI-ESM1-2-LR
## 3 3 CMCC-ESM2
## 4 4 MIROC-ES2L
## 5 5 FIO-ESM-2-0
## 6 6 CNRM-CM6-1-HR
## 7 7 INM-CM5-0
## 8 8 MRI-ESM2-0
## 9 9 EC-Earth3-Veg-LR
## 10 10 GFDL-ESM4
## 11 11 MIROC6
## 12 12 CNRM-ESM2-1
## 13 13 IPSL-CM6A-LR
## 14 14 NESM3
## 15 15 FGOALS-g3
## 16 16 ACCESS-CM2
## 17 17 AWI-CM-1-1-MR
## 18 18 CanESM5
There are no tasks for this chapter.
chill_hist_scenario_list <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_305_59")
actual_chill <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv")
chill_future_scenario_list <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_futurechill")
chills <- make_climate_scenario(
chill_hist_scenario_list,
caption = "Historic",
historic_data = actual_chill,
time_series = TRUE)
SSPs <- c("ssp126", "ssp245", "ssp370", "ssp585")
Times <- c(2050, 2085)
list_ssp <-
strsplit(names(chill_future_scenario_list), '\\.') %>%
map(2) %>%
unlist()
list_gcm <-
strsplit(names(chill_future_scenario_list), '\\.') %>%
map(3) %>%
unlist()
list_time <-
strsplit(names(chill_future_scenario_list), '\\.') %>%
map(4) %>%
unlist()
for(SSP in SSPs)
for(Time in Times)
{
# find all scenarios for the ssp and time
chill <- chill_future_scenario_list[list_ssp == SSP & list_time == Time]
names(chill) <- list_gcm[list_ssp == SSP & list_time == Time]
if(SSP == "ssp126") SSPcaption <- "SSP1"
if(SSP == "ssp245") SSPcaption <- "SSP2"
if(SSP == "ssp370") SSPcaption <- "SSP3"
if(SSP == "ssp585") SSPcaption <- "SSP5"
if(Time == "2050") Time_caption <- "2050"
if(Time == "2085") Time_caption <- "2085"
chills <- chill %>%
make_climate_scenario(
caption = c(SSPcaption,
Time_caption),
add_to = chills)
}
plot_climate_scenarios(
climate_scenario_list = chills,
metric = "Chill_Portions",
metric_label = "Chill (Chill Portions)",
texcex = 1)
## [[1]]
## [1] "time series labels"
##
## [[2]]
## code Label
## 1 1 INM-CM4-8
## 2 2 MPI-ESM1-2-LR
## 3 3 CMCC-ESM2
## 4 4 MIROC-ES2L
## 5 5 FIO-ESM-2-0
## 6 6 CNRM-CM6-1-HR
## 7 7 INM-CM5-0
## 8 8 MRI-ESM2-0
## 9 9 EC-Earth3-Veg-LR
## 10 10 GFDL-ESM4
## 11 11 MIROC6
## 12 12 CNRM-ESM2-1
## 13 13 IPSL-CM6A-LR
## 14 14 NESM3
## 15 15 FGOALS-g3
## 16 16 ACCESS-CM2
## 17 17 AWI-CM-1-1-MR
## 18 18 CanESM5
##
## [[3]]
## code Label
## 1 1 ACCESS-CM2
## 2 2 CMCC-ESM2
## 3 3 FIO-ESM-2-0
## 4 4 CNRM-CM6-1-HR
## 5 5 AWI-CM-1-1-MR
## 6 6 EC-Earth3-Veg-LR
## 7 7 GFDL-ESM4
## 8 8 CNRM-ESM2-1
## 9 9 CanESM5
## 10 10 FGOALS-g3
## 11 11 INM-CM4-8
## 12 12 INM-CM5-0
## 13 13 IPSL-CM6A-LR
## 14 14 MIROC-ES2L
## 15 15 MIROC6
## 16 16 MPI-ESM1-2-LR
## 17 17 MRI-ESM2-0
## 18 18 NESM3
##
## [[4]]
## code Label
## 1 1 FGOALS-g3
## 2 2 INM-CM4-8
## 3 3 MPI-ESM1-2-LR
## 4 4 CMCC-ESM2
## 5 5 EC-Earth3-CC
## 6 6 MIROC-ES2L
## 7 7 FIO-ESM-2-0
## 8 8 CNRM-CM6-1-HR
## 9 9 INM-CM5-0
## 10 10 MRI-ESM2-0
## 11 11 EC-Earth3-Veg-LR
## 12 12 GFDL-ESM4
## 13 13 MIROC6
## 14 14 CNRM-ESM2-1
## 15 15 IPSL-CM6A-LR
## 16 16 NESM3
## 17 17 ACCESS-CM2
## 18 18 AWI-CM-1-1-MR
##
## [[5]]
## code Label
## 1 1 FGOALS-g3
## 2 2 ACCESS-CM2
## 3 3 CMCC-ESM2
## 4 4 EC-Earth3-CC
## 5 5 FIO-ESM-2-0
## 6 6 CNRM-CM6-1-HR
## 7 7 AWI-CM-1-1-MR
## 8 8 EC-Earth3-Veg-LR
## 9 9 GFDL-ESM4
## 10 10 CNRM-ESM2-1
## 11 11 INM-CM4-8
## 12 12 INM-CM5-0
## 13 13 IPSL-CM6A-LR
## 14 14 MIROC-ES2L
## 15 15 MIROC6
## 16 16 MPI-ESM1-2-LR
## 17 17 MRI-ESM2-0
## 18 18 NESM3
##
## [[6]]
## code Label
## 1 1 CNRM-ESM2-1
## 2 2 IPSL-CM6A-LR
## 3 3 FGOALS-g3
## 4 4 EC-Earth3-AerChem
## 5 5 INM-CM4-8
## 6 6 MPI-ESM1-2-LR
## 7 7 MIROC-ES2L
## 8 8 CNRM-CM6-1-HR
## 9 9 INM-CM5-0
## 10 10 MRI-ESM2-0
## 11 11 EC-Earth3-Veg-LR
## 12 12 GFDL-ESM4
## 13 13 MIROC6
## 14 14 ACCESS-CM2
## 15 15 AWI-CM-1-1-MR
##
## [[7]]
## code Label
## 1 1 CNRM-ESM2-1
## 2 2 FGOALS-g3
## 3 3 EC-Earth3-AerChem
## 4 4 ACCESS-CM2
## 5 5 CNRM-CM6-1-HR
## 6 6 AWI-CM-1-1-MR
## 7 7 EC-Earth3-Veg-LR
## 8 8 GFDL-ESM4
## 9 9 INM-CM4-8
## 10 10 INM-CM5-0
## 11 11 IPSL-CM6A-LR
## 12 12 MIROC-ES2L
## 13 13 MIROC6
## 14 14 MPI-ESM1-2-LR
## 15 15 MRI-ESM2-0
##
## [[8]]
## code Label
## 1 1 CIESM
## 2 2 NESM3
## 3 3 CNRM-ESM2-1
## 4 4 IPSL-CM6A-LR
## 5 5 FGOALS-g3
## 6 6 CMCC-ESM2
## 7 7 INM-CM4-8
## 8 8 MPI-ESM1-2-LR
## 9 9 EC-Earth3-CC
## 10 10 FIO-ESM-2-0
## 11 11 MIROC-ES2L
## 12 12 CNRM-CM6-1-HR
## 13 13 INM-CM5-0
## 14 14 MRI-ESM2-0
## 15 15 EC-Earth3-Veg-LR
## 16 16 GFDL-ESM4
## 17 17 MIROC6
## 18 18 ACCESS-CM2
## 19 19 AWI-CM-1-1-MR
##
## [[9]]
## code Label
## 1 1 CIESM
## 2 2 CNRM-ESM2-1
## 3 3 FGOALS-g3
## 4 4 CMCC-ESM2
## 5 5 ACCESS-CM2
## 6 6 EC-Earth3-CC
## 7 7 FIO-ESM-2-0
## 8 8 CNRM-CM6-1-HR
## 9 9 AWI-CM-1-1-MR
## 10 10 EC-Earth3-Veg-LR
## 11 11 GFDL-ESM4
## 12 12 INM-CM4-8
## 13 13 INM-CM5-0
## 14 14 IPSL-CM6A-LR
## 15 15 MIROC-ES2L
## 16 16 MIROC6
## 17 17 MPI-ESM1-2-LR
## 18 18 MRI-ESM2-0
## 19 19 NESM3
for(nam in names(chills[[1]]$data))
{
ch <- chills[[1]]$data[[nam]]
ch[,"GCM"] <- "none"
ch[,"SSP"] <- "none"
ch[,"Year"] <- as.numeric(nam)
if(nam == names(chills[[1]]$data)[1])
past_simulated <- ch else
past_simulated <- rbind(past_simulated,
ch)
}
past_simulated["Scenario"] <- "Historical"
head(past_simulated)
## Season End_year Season_days Data_days Perc_complete Chill_Portions
## 1 2001/2002 2002 120 120 100 85.93502
## 2 2002/2003 2003 120 120 100 84.11282
## 3 2003/2004 2004 120 120 100 83.33690
## 4 2004/2005 2005 121 121 100 88.70952
## 5 2005/2006 2006 120 120 100 89.60029
## 6 2006/2007 2007 120 120 100 84.94518
## GDH Frost_H GCM SSP Year Scenario
## 1 3996.664 450 none none 1990 Historical
## 2 3344.907 479 none none 1990 Historical
## 3 1889.238 721 none none 1990 Historical
## 4 2464.805 361 none none 1990 Historical
## 5 3750.400 214 none none 1990 Historical
## 6 3091.575 432 none none 1990 Historical
past_simulated <- past_simulated %>% filter(Perc_complete == 100)
past_observed <- chills[[1]][["historic_data"]]
head(past_observed)
## X Season End_year Season_days Data_days Perc_complete Chill_Portions
## 1 1 1990/1991 1991 120 120 100 77.46451
## 2 2 1991/1992 1992 120 120 100 84.12147
## 3 3 1992/1993 1993 121 121 100 85.89337
## 4 4 1993/1994 1994 120 120 100 83.83349
## 5 5 1994/1995 1995 120 120 100 86.00229
## 6 6 1995/1996 1996 120 120 100 67.05746
## GDH Frost_H
## 1 1964.757 839
## 2 2481.648 491
## 3 3143.898 482
## 4 1904.884 591
## 5 5425.112 317
## 6 1612.944 1224
for(i in 2:length(chills))
for(nam in names(chills[[i]]$data))
{ch <- chills[[i]]$data[[nam]]
ch[,"GCM"] <- nam
ch[,"SSP"] <- chills[[i]]$caption[1]
ch[,"Year"] <- chills[[i]]$caption[2]
if(i == 2 & nam == names(chills[[i]]$data)[1])
future_data <- ch else
future_data <- rbind(future_data,ch)
}
head(future_data)
## Season End_year Season_days Data_days Perc_complete Chill_Portions
## 1 2001/2002 2002 120 120 100 86.99319
## 2 2002/2003 2003 120 120 100 86.87197
## 3 2003/2004 2004 120 120 100 85.62412
## 4 2004/2005 2005 121 121 100 90.50165
## 5 2005/2006 2006 120 120 100 90.64776
## 6 2006/2007 2007 120 120 100 86.32296
## GDH Frost_H GCM SSP Year
## 1 5318.501 352 INM-CM4-8 SSP1 2050
## 2 4400.557 409 INM-CM4-8 SSP1 2050
## 3 2744.363 578 INM-CM4-8 SSP1 2050
## 4 3436.639 242 INM-CM4-8 SSP1 2050
## 5 5051.034 111 INM-CM4-8 SSP1 2050
## 6 4095.937 350 INM-CM4-8 SSP1 2050
metric <- "GDH"
axis_label <- "Heat (in GDH)"
rng <- range(past_observed[[metric]],
past_simulated[[metric]],
future_data[[metric]])
rng
## [1] 1088.303 19014.867
past_plot <- ggplot() +
geom_boxplot(data = past_simulated,
aes_string("as.numeric(Year)",
metric,
group = "Year"),
fill = "skyblue")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
past_plot <- past_plot +
scale_y_continuous(
limits = c(0,
round(rng[2] + rng[2]/10))) +
labs(x = "Year",
y = axis_label)
past_plot <- past_plot +
facet_grid(~ Scenario) +
theme_bw(base_size = 15)
past_plot <- past_plot +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 45,
hjust = 1))
# add historic data
past_plot <- past_plot +
geom_point(data = past_observed,
aes_string("End_year",
metric),
col = "blue")
past_plot
y <- 2050
future_2050 <-
ggplot(data = future_data[future_data$Year == y,]) +
geom_boxplot(aes_string("GCM",
metric,
fill = "GCM"))
future_2050 <- future_2050 +
facet_wrap(vars(SSP), nrow = 1) +
scale_x_discrete(labels = NULL,
expand = expansion(add = 1))
library(ggpmisc)
future_2050 <- future_2050 +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center",
npcy = "top",
label = Year),
size = 5)
future_2050 <- future_2050 +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0,
0,
0,
0,
"cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
future_2050
future_plot_list <- list()
time_points <- c(2050, 2085)
for(y in time_points)
{
future_plot_list[[which(y == time_points)]] <-
ggplot(data = future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM",
metric,
fill="GCM")) +
facet_wrap(vars(SSP), nrow = 1) +
scale_x_discrete(labels = NULL,
expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center",
npcy = "top",
label = Year),
size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0,
0,
0,
0,
"cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(
hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
}
library(patchwork)
both_plots <- past_plot + future_plot_list
plot <- both_plots +
plot_layout(guides = "collect",
widths = c(1,
rep(1.8,
length(future_plot_list))))
plot <- plot & theme(legend.position = "bottom",
legend.text = element_text(size = 8),
legend.title = element_text(size = 10),
axis.title.x = element_blank())
plot
metric <- "Chill_Portions"
axis_label <- "Chill (in CP)"
rng <- range(past_observed[[metric]],
past_simulated[[metric]],
future_data[[metric]])
past_plot <- ggplot() +
geom_boxplot(data = past_simulated,
aes_string("as.numeric(Year)",
metric,group="Year"),
fill="skyblue") +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
labs(x = "Year", y = axis_label) +
facet_grid(~ Scenario) +
theme_bw(base_size = 15) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45, hjust=1)) +
geom_point(data = past_observed,
aes_string("End_year",
metric),
col="blue")
future_plot_list <- list()
for(y in c(2050,
2085))
{
future_plot_list[[which(y == c(2050,2085))]] <-
ggplot(data = future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM",
metric,
fill="GCM")) +
facet_wrap(vars(SSP), nrow = 1) +
scale_x_discrete(labels = NULL,
expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center",
npcy = "top",
label = Year),
size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0,
0,
0,
0,
"cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
}
plot <- (past_plot +
future_plot_list +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
) & theme(legend.position = "bottom",
legend.text = element_text(size = 8),
legend.title = element_text(size = 10),
axis.title.x=element_blank())
plot
metric <- "Frost_H"
axis_label <- "Frost duration (in hours)"
# get extreme values for the axis scale
rng <- range(past_observed[[metric]],
past_simulated[[metric]],
future_data[[metric]])
past_plot <- ggplot() +
geom_boxplot(data = past_simulated,
aes_string("as.numeric(Year)",
metric,group="Year"),
fill="skyblue") +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
labs(x = "Year", y = axis_label) +
facet_grid(~ Scenario) +
theme_bw(base_size = 15) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45, hjust=1)) +
geom_point(data = past_observed,
aes_string("End_year",
metric),
col="blue")
future_plot_list <- list()
for(y in c(2050,
2085))
{
future_plot_list[[which(y == c(2050,2085))]] <-
ggplot(data = future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM",
metric,
fill="GCM")) +
facet_wrap(vars(SSP), nrow = 1) +
scale_x_discrete(labels = NULL,
expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center",
npcy = "top",
label = Year),
size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0,
0,
0,
0,
"cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
}
plot <- (past_plot +
future_plot_list +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
) & theme(legend.position = "bottom",
legend.text = element_text(size = 8),
legend.title = element_text(size = 10),
axis.title.x=element_blank())
plot
plot_scenarios_gg <- function(past_observed,
past_simulated,
future_data,
metric,
axis_label,
time_points)
{
rng <- range(past_observed[[metric]],
past_simulated[[metric]],
future_data[[metric]])
past_plot <- ggplot() +
geom_boxplot(data = past_simulated,
aes_string("as.numeric(Year)",
metric,
group="Year"),
fill="skyblue") +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
labs(x = "Year", y = axis_label) +
facet_grid(~ Scenario) +
theme_bw(base_size = 15) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45,
hjust=1)) +
geom_point(data = past_observed,
aes_string("End_year",
metric),
col="blue")
future_plot_list <- list()
for(y in time_points)
{
future_plot_list[[which(y == time_points)]] <-
ggplot(data = future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM",
metric,
fill="GCM")) +
facet_wrap(vars(SSP), nrow = 1) +
scale_x_discrete(labels = NULL,
expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0,
round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center",
npcy = "top",
label = Year),
size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0,
0,
0,
0,
"cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
}
plot <- (past_plot +
future_plot_list +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
) & theme(legend.position = "bottom",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
}
plot_scenarios_gg(past_observed = past_observed,
past_simulated = past_simulated,
future_data = future_data,
metric = "GDH",
axis_label = "Heat (in Growing Degree Hours)",
time_points = c(2050, 2085))
ggsave("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_heat_plot.png",
width = 10,
height = 8,
dpi = 600)
plot_scenarios_gg(past_observed = past_observed,
past_simulated = past_simulated,
future_data = future_data,
metric = "Chill_Portions",
axis_label = "Chill (in Chill Portions)",
time_points = c(2050, 2085))
ggsave("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_Chill_portions_plot.png",
width = 10,
height = 8,
dpi = 600)
plot_scenarios_gg(past_observed = past_observed,
past_simulated = past_simulated,
future_data = future_data,
metric = "Frost_H",
axis_label = "Frost duration (in hours)",
time_points = c(2050, 2085))
ggsave("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_frost_plot.png",
width = 10,
height = 8,
dpi = 600)
hourly_models <- list(Chilling_units = chilling_units,
Low_chill = low_chill_model,
Modified_Utah = modified_utah_model,
North_Carolina = north_carolina_model,
Positive_Utah = positive_utah_model,
Chilling_Hours = Chilling_Hours,
Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model)
daily_models <- list(Rate_of_Chill = rate_of_chill,
Chill_Days = chill_days,
Exponential_Chill = exponential_chill,
# Triangular_Chill_Haninnen = triangular_chill_1,
Triangular_Chill_Legave = triangular_chill_2)
metrics <- c(names(daily_models),
names(hourly_models))
model_labels = c("Rate of Chill",
"Chill Days",
"Exponential Chill",
# "Triangular Chill (Häninnen)",
"Triangular Chill (Legave)",
"Chilling Units",
"Low-Chill Chill Units",
"Modified Utah Chill Units",
"North Carolina Chill Units",
"Positive Utah Chill Units",
"Chilling Hours",
"Utah Chill Units",
"Chill Portions")
data.frame(Metric = model_labels, 'Function name' = metrics)
## Metric Function.name
## 1 Rate of Chill Rate_of_Chill
## 2 Chill Days Chill_Days
## 3 Exponential Chill Exponential_Chill
## 4 Triangular Chill (Legave) Triangular_Chill_Legave
## 5 Chilling Units Chilling_units
## 6 Low-Chill Chill Units Low_chill
## 7 Modified Utah Chill Units Modified_Utah
## 8 North Carolina Chill Units North_Carolina
## 9 Positive Utah Chill Units Positive_Utah
## 10 Chilling Hours Chilling_Hours
## 11 Utah Chill Units Utah_Chill_Units
## 12 Chill Portions Chill_Portions
Wesel_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")
Temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_scenarios")
Start_JDay <- 305
End_JDay <- 59
daily_models_past_scenarios <-
tempResponse_list_daily(Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_past_scenarios <- lapply(
daily_models_past_scenarios,
function(x) x[which(x$Perc_complete>90),])
hourly_models_past_scenarios<-
tempResponse_daily_list(Temps,
latitude = 51.405,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models = hourly_models,
misstolerance = 10)
past_scenarios <- daily_models_past_scenarios
past_scenarios <- lapply(
names(past_scenarios),
function(x)
cbind(past_scenarios[[x]],
hourly_models_past_scenarios[[x]][,names(hourly_models)]))
names(past_scenarios) <- names(daily_models_past_scenarios)
daily_models_observed <-
tempResponse_daily(Wesel_temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models = daily_models)
daily_models_observed <-
daily_models_observed[which(daily_models_observed$Perc_complete>90),]
hourly_models_observed <-
tempResponse_daily_list(Wesel_temps,
latitude=51.405,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models = hourly_models,
misstolerance = 10)
past_observed <- cbind(
daily_models_observed,
hourly_models_observed[[1]][,names(hourly_models)])
save_temperature_scenarios(past_scenarios,
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_multichill_305_59_historic")
write.csv(past_observed,
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_multichill_305_59_observed.csv",
row.names=FALSE)
SSPs <- c("ssp126", "ssp245","ssp370", "ssp585")
Times <- c(2050, 2085)
future_temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate","Wesel_future_")
list_ssp <-
strsplit(names(future_temps), '\\.') %>%
map(2) %>%
unlist()
list_gcm <-
strsplit(names(future_temps), '\\.') %>%
map(3) %>%
unlist()
list_time <-
strsplit(names(future_temps), '\\.') %>%
map(4) %>%
unlist()
for(SSP in SSPs)
{
for(Time in Times)
{
Temps <- future_temps[list_ssp == SSP & list_time == Time]
names(Temps) <- list_gcm[list_ssp == SSP & list_time == Time]
daily_models_future_scenarios <- tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models = daily_models)
daily_models_future_scenarios<-lapply(
daily_models_future_scenarios,
function(x) x[which(x$Perc_complete>90),])
hourly_models_future_scenarios<-
tempResponse_daily_list(
Temps,
latitude = 51.405,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
future_scenarios <- daily_models_future_scenarios
future_scenarios <- lapply(
names(future_scenarios),
function(x)
cbind(future_scenarios[[x]],
hourly_models_future_scenarios[[x]][,names(hourly_models)]))
names(future_scenarios)<-names(daily_models_future_scenarios)
chill<-future_scenarios
save_temperature_scenarios(
chill,
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
paste0("Wesel_multichill_305_59_",Time,"_",SSP))
}
}
chill_past_scenarios <- load_temperature_scenarios(
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
"Wesel_multichill_305_59_historic")
chill_observed <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_multichill_305_59_observed.csv")
chills <- make_climate_scenario(chill_past_scenarios,
caption = "Historical",
historic_data = chill_observed,
time_series = TRUE)
for(SSP in SSPs)
for(Time in Times)
{
chill <- load_temperature_scenarios(
"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
paste0("Wesel_multichill_305_59_",Time,"_",SSP))
if(SSP == "ssp126") SSPcaption <- "SSP1"
if(SSP == "ssp245") SSPcaption <- "SSP2"
if(SSP == "ssp370") SSPcaption <- "SSP3"
if(SSP == "ssp585") SSPcaption <- "SSP5"
if(Time == "2050") Time_caption <- "2050"
if(Time == "2085") Time_caption <- "2085"
chills <- make_climate_scenario(chill,
caption = c(SSPcaption,
Time_caption),
add_to = chills)
}
for(i in 1:length(chills))
{ch <- chills[[i]]
if(ch$caption[1] == "Historical")
{GCMs <- rep("none",length(names(ch$data)))
SSPs <- rep("none",length(names(ch$data)))
Years <- as.numeric(ch$labels)
Scenario <- rep("Historical",
length(names(ch$data)))} else
{GCMs <- names(ch$data)
SSPs <- rep(ch$caption[1],
length(names(ch$data)))
Years <- rep(as.numeric(ch$caption[2]),
length(names(ch$data)))
Scenario <- rep("Future",
length(names(ch$data)))}
for(nam in names(ch$data))
{for(met in metrics)
{temp_res <-
data.frame(Metric = met,
GCM = GCMs[which(nam == names(ch$data))],
SSP = SSPs[which(nam == names(ch$data))],
Year = Years[which(nam == names(ch$data))],
Result = quantile(ch$data[[nam]][,met],0.1),
Scenario = Scenario[which(nam == names(ch$data))])
if(i == 1 & nam == names(ch$data)[1] & met == metrics[1])
results <- temp_res else
results <- rbind(results,
temp_res)
}
}
}
for(met in metrics)
results[which(results$Metric == met),"SWC"] <-
results[which(results$Metric == met),"Result"]/
results[which(results$Metric == met & results$Year == 1990),
"Result"]-1
rng = range(results$SWC)
p_future <- ggplot(results[which(!results$GCM == "none"),],
aes(GCM,
y = factor(Metric,
levels = metrics),
fill = SWC)) +
geom_tile()
p_future <-
p_future +
facet_grid(SSP ~ Year)
p_future <-
p_future +
theme_bw(base_size = 15) +
theme(axis.text = element_text(size=6))
library(colorRamps)
p_future <-
p_future +
scale_fill_gradientn(colours = matlab.like(15),
labels = scales::percent,
limits = rng)
p_future <-
p_future +
theme(axis.text.x = element_text(angle = 75,
hjust = 1,
vjust = 1)) +
labs(fill = "Change in\nSafe Winter Chill\nsince 1980") +
scale_y_discrete(labels = model_labels) +
ylab("Chill metric")
p_future
p_past<-
ggplot(results[which(results$GCM == "none"),],
aes(Year,
y = factor(Metric,
levels=metrics),
fill = SWC)) +
geom_tile()
p_past<-
p_past +
theme_bw(base_size = 15) +
theme(axis.text = element_text(size = 6))
p_past<-
p_past +
scale_fill_gradientn(colours = matlab.like(15),
labels = scales::percent,
limits = rng)
p_past<-
p_past +
scale_x_continuous(position = "top")
p_past<-
p_past +
labs(fill = "Change in\nSafe Winter Chill\nsince 1980") +
scale_y_discrete(labels = model_labels) +
ylab("Chill metric")
p_past
chill_comp_plot<-
(p_past +
p_future +
plot_layout(guides = "collect",
nrow = 2,
heights = c(1,3))) &
theme(legend.position = "right",
strip.background = element_blank(),
strip.text = element_text(face = "bold"))
chill_comp_plot
hist_results <- results[which(results$GCM == "none"),]
hist_results$SSP <- "SSP1"
hist_results_2 <- hist_results
hist_results_2$SSP <- "SSP2"
hist_results_3 <- hist_results
hist_results_3$SSP <- "SSP3"
hist_results_4 <- hist_results
hist_results_4$SSP <- "SSP5"
hist_results <- rbind(hist_results,
hist_results_2,
hist_results_3,
hist_results_4)
future_results <- results[which(!results$GCM == "none"),]
GCM_aggregate <- aggregate(
future_results$SWC,
by=list(future_results$Metric,
future_results$SSP,
future_results$Year),
FUN=mean)
colnames(GCM_aggregate) <- c("Metric",
"SSP",
"Year",
"SWC")
SSP_Time_series<-rbind(hist_results[,c("Metric",
"SSP",
"Year",
"SWC")],
GCM_aggregate)
SSP_Time_series$Year <- as.numeric(SSP_Time_series$Year)
chill_change_plot<-
ggplot(data = SSP_Time_series,
aes(x = Year,
y = SWC,
col = factor(Metric,
levels = metrics))) +
geom_line(lwd = 1.3) +
facet_wrap(~SSP,
nrow = 4) +
theme_bw(base_size = 10) +
labs(col = "Change in\nSafe Winter Chill\nsince 1980") +
scale_color_discrete(labels = model_labels) +
scale_y_continuous(labels = scales::percent) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold")) +
ylab("Safe Winter Chill")
chill_change_plot
library(gganimate)
library(gifski)
## Warning: Paket 'gifski' wurde unter R Version 4.4.2 erstellt
library(png)
library(transformr)
## Warning: Paket 'transformr' wurde unter R Version 4.4.2 erstellt
ccp<-chill_change_plot +
transition_reveal(Year)
animate(ccp, fps = 10)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
anim_save("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/chill_comparison_animation.gif",
animation = last_animation())